//
// Created by 王岩松 on 2020/3/8.
//

#ifndef LPC_MATRIX_H
#define LPC_MATRIX_H

#include <vector>
#include <cassert>
#include <ostream>
#include <cmath>

class matrix{
private:
    double EPS = 1e-10;
    bool sign(double a) {
        return a < -EPS ? -1 : (a > EPS);
    }
    void checkZero() {
        for(int i = 0; i < n; ++ i) {
            for(int j = 0; j < m; ++ j) {
                if(sign(A[i][j])==0) A[i][j] = 0;
            }
        }
    }
public:
    int n,m;
    std::vector<std::vector<double>> A;
    matrix();
    matrix(int _n,int _m);
    matrix(int _n,int _m,double val);
    matrix(std::vector<std::vector<double>> _A);
    matrix operator + (matrix a) const {
        assert(a.n == n && a.m == m);
        matrix ans(n,m);
        for(int i = 0; i < n; ++ i) {
            for(int j = 0;  j < m; ++ j) {
                ans.A[i][j]  = A[i][j] + a.A[i][j];
            }
        }
        return ans;
    }
    matrix operator - (matrix a)  const {
        assert(a.n == n && a.m == m);
        matrix ans(n,m);
        for(int i = 0; i < n; ++ i) {
            for(int j = 0;  j < m; ++ j) {
                ans.A[i][j]  = A[i][j] - a.A[i][j];
            }
        }
        return ans;
    }
    matrix operator * (matrix a) const {
        assert(m == a.n);
        matrix ans(n,a.m);
        for(int k = 0; k < m; ++ k) {
            for(int i = 0; i < n; ++ i) {
                for(int j = 0;  j < a.m; ++ j) {
                    ans.A[i][j] += A[i][k] * a.A[k][j];
                }
            }
        }
        return ans;
    }
    matrix operator / (matrix a)  const{
        return a.inverse()*(*this);
    }

    matrix operator + (double a) const {
        matrix ans(n,m);
        for(int i = 0; i < n; ++ i) {
            for(int j = 0; j < m; ++ j) {
                ans[i][j] = A[i][j] + a;
            }
        }
        return ans;
    }
    matrix operator - (double a) const {
        matrix ans(n,m);
        for(int i = 0; i < n; ++ i) {
            for(int j = 0; j < m; ++ j) {
                ans[i][j] = A[i][j] - a;
            }
        }
        return ans;
    }
    matrix operator * (double a) const {
        matrix ans(n,m);
        for(int i = 0; i < n; ++ i) {
            for(int j = 0; j < m; ++ j) {
                ans[i][j] = A[i][j] * a;
            }
        }
        return ans;
    }
    matrix operator / (double a) const {
        assert(a != 0);
        matrix ans(n,m);
        for(int i = 0; i < n; ++ i) {
            for(int j = 0; j < m; ++ j) {
                ans[i][j] = A[i][j] / a;
            }
        }
        return ans;
    }
    std::vector<double>& operator [](int index) {
        return A[index];
    }
    matrix dot(matrix a){
        assert(n == a.n && m == a.m);
        matrix ans(n,m);
        for(int i = 0; i < n; ++ i) {
            for(int j = 0; j < m; ++ j) {
                ans[i][j] = A[i][j] + a[i][j];
            }
        }
        return ans;
    }
    matrix appendRight(matrix a);   //在右边扩展一个矩阵
    matrix appendDown(matrix a);    //在下边扩展一个矩阵
    matrix transpose();             //矩阵的转置
    matrix gauss();                 //对角线gauss消元成1
    matrix inverse();               //矩阵的逆
    void Free();
};

#endif //LPC_MATRIX_H
